Simulation-efficient shortest probability intervals

Bayesian highest posterior density (HPD) intervals can be estimated directly from simultions via empirical shortest intervals. Unfortunately, these can be noisy (that is, have a high Monte Carlo error). In this paper, we propose an optimal weighting strategy using quadratic programming to obtain a more computationally stable HPD, or in general, shortest probability interval (Spin). We prove the consistency of our method. Simulation studies show that Spin has lower Monte Carlo error than empirical shortest intervals. We also apply our approach to two applied Bayesian examples featuring highlyasymmetric posterior distributions. We implement the new method in an R package (SPIn) so it can be routinely used in post-processing of Bayesian simulations.